In [1]:
from util import *
from glob import glob

/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:124: UserWarning: The Shapely GEOS version (3.11.2-CAPI-1.17.2) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow.
  warnings.warn(
In [2]:
df = load_AOIs()
df
Out[2]:
Taranaki AOI SSP 4.5 (p50) SSP 4.5 (p83) SSP 8.5 (p50) SSP 8.5 (p83) Rate SSP 4.5 (p50) Rate SSP 4.5 (p83) Rate SSP 8.5 (p50) Rate SSP 8.5 (p83) match match_score
7 NORTH TongaporutuRiver 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 TongapurutuRiverCliffs 93.750000
11 SOUTH HangatahuaRiver_South 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 HangatahuRiver_South 97.435897
21 SOUTH Rahotu 0.58 0.78 0.84 1.10 0.0058 0.0078 0.0084 0.0110 Rahotu 100.000000
20 SOUTH Pihama 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Pihama 100.000000
19 SOUTH OpunakeBeach 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OpunakeBeachCliffs 100.000000
18 SOUTH OhaweBeach 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 OhaweBeach 100.000000
17 SOUTH Oeo 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Oeo 100.000000
16 SOUTH Manutahi 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Manutahi 100.000000
15 SOUTH ManaBay 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 ManaBayCliffs 100.000000
14 SOUTH KaupokonuiBeach 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 KaupokonuiBeach 100.000000
13 SOUTH Kakaramea 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Kakaramea 100.000000
12 SOUTH Hawera_WaihiBeach 0.57 0.78 0.83 1.10 0.0057 0.0078 0.0083 0.0110 Hawera_WaihiBeach 100.000000
0 NORTH Mohakatino_River 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 MohakatinoRiver 100.000000
10 SOUTH CapeEgmont 0.58 0.78 0.84 1.10 0.0058 0.0078 0.0084 0.0110 CapeEgmont 100.000000
9 NORTH Waitara 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Waitara 100.000000
8 NORTH UrenuiRiver_North 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 UrenuiRiverNorthCliffs 100.000000
6 NORTH PariokariwaPoint 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 PariokariwaPointCliffs 100.000000
5 NORTH Onaero 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OnaeroCliff 100.000000
4 NORTH Oakura_South 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 OakuraSouth 100.000000
3 NORTH Oakura 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 Oakura 100.000000
2 NORTH NewPlymouth_Airport 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 NewPlymouthAirport 100.000000
1 NORTH NewPlymouth 0.57 0.78 0.84 1.10 0.0057 0.0078 0.0084 0.0110 NewPlymouthCliffs 100.000000
22 SOUTH WainuiBeach 0.57 0.78 0.83 1.09 0.0057 0.0078 0.0083 0.0109 WainuiBeach 100.000000
23 SOUTH WaipipiBeach 0.57 0.78 0.83 1.09 0.0057 0.0078 0.0083 0.0109 WaipipiBeachCliffs 100.000000
In [3]:
site = df.match.sample(1).iloc[0]
site
Out[3]:
'Oeo'
In [4]:
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
if site == "ManaBayCliffs":
  print("Flipping")
  for k, v in transect_metadata.items():
    transect_metadata[k]["Azimuth"] = v["Azimuth"] + 180

/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
In [5]:
linear_models = fit(gdf, transect_metadata)
linear_models.loc[linear_models.slope > 0, "slope"] = pd.NA
rolled_slopes = linear_models.groupby("group").slope.rolling(10, min_periods=1).mean().dropna().reset_index(level=0)
linear_models.slope = rolled_slopes.slope
linear_models.dropna(inplace=True)
linear_models
Out[5]:
TransectID slope intercept group
0 2 -0.025895 11.839523 0
1 3 -0.026103 11.394666 0
2 4 -0.032111 12.067457 0
3 5 -0.037306 12.952147 0
4 6 -0.041094 14.552892 0
... ... ... ... ...
665 806 -0.088458 17.144344 50
666 807 -0.098385 17.888823 50
667 808 -0.100127 16.604210 50
668 809 -0.105771 5.736368 50
669 810 -0.105771 4.466733 50

657 rows × 4 columns

In [6]:
results = predict(gdf, linear_models, transect_metadata)
results
Out[6]:
TransectID BaselineID group Year ocean_point linear_model_point sqrt_model_point BH_model_point Sunamura_model_point
0 2.0 1 0.0 2100 POINT (1687908.5767938623 5620307.34597674) POINT (1688034.2711515864 5620793.369995707) POINT (1688035.5534006841 5620798.3280850835) POINT (1688043.1672309272 5620827.768582187) POINT (1688037.0176938204 5620803.9900865285)
1 3.0 1 0.0 2100 POINT (1687915.1713428975 5620306.442715614) POINT (1688044.0861092217 5620791.621524045) POINT (1688045.4118067573 5620796.610869867) POINT (1688053.2101313008 5620825.9603533875) POINT (1688046.9030184844 5620802.223138728)
2 4.0 1 0.0 2100 POINT (1687923.1177022296 5620304.70048816) POINT (1688053.6788356411 5620788.917562359) POINT (1688055.3321562614 5620795.049295843) POINT (1688062.928654537 5620823.222719841) POINT (1688056.5343737882 5620799.508006012)
3 5.0 1 0.0 2100 POINT (1687929.0214491673 5620302.536182147) POINT (1688063.1581929883 5620785.8900468) POINT (1688065.1311122181 5620792.999345171) POINT (1688072.6592310544 5620820.126478494) POINT (1688066.0911000874 5620796.458604978)
4 6.0 1 0.0 2100 POINT (1687943.6634779107 5620297.272404213) POINT (1688072.646862405 5620782.368306764) POINT (1688074.7352682156 5620790.222629285) POINT (1688081.7768776012 5620816.705543141) POINT (1688075.4651071057 5620792.967493336)
... ... ... ... ... ... ... ... ... ...
652 806.0 1 50.0 2100 POINT (1687794.0844616203 5620310.788355438) POINT (1687889.1002257578 5620803.734863906) POINT (1687892.4113368606 5620820.913072934) POINT (1687895.824934831 5620838.622988925) POINT (1687891.1748046514 5620814.497882353)
653 807.0 1 50.0 2100 POINT (1687821.543196849 5620306.082154956) POINT (1687899.1019846308 5620802.495838878) POINT (1687902.105572404 5620821.720250634) POINT (1687904.5866393007 5620837.600276526) POINT (1687900.7938008038 5620813.324279123)
654 808.0 1 50.0 2100 POINT (1687873.6267297296 5620298.320708633) POINT (1687909.2520916292 5620802.726387996) POINT (1687910.64721948 5620822.479462039) POINT (1687911.7553046227 5620838.168409724) POINT (1687910.0242252424 5620813.6587283)
655 809.0 1 50.0 2100 POINT (1687893.8018985782 5620296.623321697) POINT (1687919.8082378658 5620813.567189055) POINT (1687920.859272824 5620834.459251945) POINT (1687921.5934341024 5620849.052623445) POINT (1687920.3588562713 5620824.512166921)
656 810.0 1 50.0 2100 POINT (1687939.3591070531 5620298.028960121) POINT (1687928.690593155 5620815.9733632095) POINT (1687928.2598098365 5620836.8874109285) POINT (1687927.9589022126 5620851.496139112) POINT (1687928.4649134837 5620826.929858534)

657 rows × 9 columns

In [7]:
gpd.GeoSeries(results.linear_model_point).distance(gpd.GeoSeries(results.Sunamura_model_point)).describe()
Out[7]:
count    657.000000
mean      10.961338
std        0.009703
min       10.924303
25%       10.957830
50%       10.964826
75%       10.967859
max       10.972655
dtype: float64
In [8]:
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
poly = prediction_results_to_polygon(results)
In [9]:
poly.explore(tiles="Esri.WorldImagery")
Out[9]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [10]:
site = df.match.sample(1).iloc[0]
site = "CapeEgmont"
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
m = gpd.read_file(f"Shapefiles/{site}_TransectLines.shp").explore("TransectID", tiles="Esri.WorldImagery", max_zoom=22)
linear_models = fit(gdf, transect_metadata)
linear_models.loc[linear_models.slope > 0, "slope"] = pd.NA
rolled_slopes = linear_models.groupby("group").slope.rolling(10, min_periods=1).mean().dropna().reset_index(level=0)
linear_models.slope = rolled_slopes.slope
linear_models.dropna(inplace=True)
results = predict(gdf, linear_models, transect_metadata)
results["geometry"] = results["linear_model_point"]
pd.concat([gdf, results])[["Year", "geometry"]].explore("Year", m=m)
gpd.GeoSeries(results.ocean_point, crs=2193).explore(m=m, color="blue")
#gpd.GeoSeries(poly, crs=2193).explore(color="pink", m=m)
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
polygon = prediction_results_to_polygon(results)
polygon.explore(m=m)
m

CapeEgmont

/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
<ipython-input-10-8352519b00ac>:15: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
  results["geometry"] = results["linear_model_point"]
Out[10]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [11]:
run_all_parallel()

/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)

Flipping

/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
In [12]:
def read_file(f):
  df = gpd.read_file(f)
  df["filename"] = f
  return df
samples = pd.concat(read_file(f) for f in glob("Projected_Shoreline_Polygons/*.shp"))
samples.explore("filename", tiles="Esri.WorldImagery", style_kwds={"fill": False})
Out[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [13]:
# Compare different SLR rates
site = df.match.sample(1).iloc[0]
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
linear_models = fit(gdf, transect_metadata)
# Erosion only
linear_models.loc[linear_models.slope > 0, "slope"] = pd.NA
rolled_slopes = linear_models.groupby("group").slope.rolling(10, min_periods=1).mean().dropna().reset_index(level=0)
linear_models.slope = rolled_slopes.slope
linear_models.dropna(inplace=True)
low_proj_slr = predict(gdf, linear_models, transect_metadata, Proj_SLR=0.007)
high_proj_slr = predict(gdf, linear_models, transect_metadata, Proj_SLR=0.015)
m = gpd.GeoSeries(low_proj_slr.sqrt_model_point, crs=2193).explore(color="blue", tiles="Esri.WorldImagery")
gpd.GeoSeries(high_proj_slr.sqrt_model_point, crs=2193).explore(color="green", m=m)

Oakura

/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
  return lib.distance(a, b, **kwargs)
Out[13]:
Make this Notebook Trusted to load map: File -> Trust Notebook